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Although species longevity is subject to a diverse range of selective forces, the mortality curves 
of a wide variety of organisms are rather similar. We argue that aging and its universal character- 
istics may have evolved by means of a gradual increase in the systemic interdependence between 
a large collection of biochemical or mechanical components. Modeling the organism as a depen- 
dency network which we create using a constructive evolutionary process, we age it by allowing 
nodes to be broken or repaired according to a probabilistic algorithm that accounts for random 
failures/repairs and dependencies. Our simulations show that the network slowly accumulates dam- 
age and then catastrophically collapses. We use our simulations to fit experimental data for the 
time dependent mortality rates of a variety of multicellular organisms and even complex machines 
such as automobiles. Our study suggests that aging is an emergent finite-size effect in networks 
with dynamical dependencies and that the qualitative and quantitative features of aging are not 
sensitively dependent on the details of system structure. 



I. INTRODUCTION 

For a collection of s radioactive atoms, the probability 
of decay is a constant, so that the fraction of atoms that 
decays per unit time —(ds/dt)/s = fi does not change in 
time. In other words, an "old" atom is equally likely to 
decay as a "young" one. Contrastingly, in complex struc- 
tures such as organizations, organisms and machines, one 
finds that the relative fraction that dies per unit time, 
the mortality or aging parameter fi(t) varies, and typi- 
cally increases in time. In living systems fi{t) increases 
exponentially (commonly known as the Gompertz Law) 
up until a "late-life plateau", after which aging deceler- 
ates. Moreover, the functional form of fi(t) for a wide 
variety of organisms is remarkably similar [IH3]. 

The origins of biological aging has been sought in 
two broad classes of theories [H [5J , which may be non- 
exclusive. The first, mechanistic approach, aims to un- 
derstand aging in terms of mechanical and biochemical 
processes such as telomere shortening [5] or reactive oxy- 
gen species damage [7]. The second approach, considers 
aging as the outcome of evolutionary forces [5] . The early 
evolutionary theories are based on the observation that 
selective pressure is larger for traits that appear earlier 
in life [9HTT]. As a result, aging, it has been argued, 
could be due to late-acting deleterious mutations accu- 
mulated over generations (mutation accumulation the- 
ory, MA) H2] or due to mutations that increase fitness 
early in life at the cost of decreasing fitness later in life 
(antagonistic pleiotropy theory, AP) [TH] |T3J 03] . Physi- 
ological variants of AP consider the relative energy cost 
of avoiding aging damage versus reproducing |151 I16j : 
Mutations diverting energy away from repair and main- 
tenance activities to earlier sexual development and high 
reproduction rate can be favored. 

A more recent collection of evolutionary theories view 
aging as a group selection effect (GS); e.g. it has been 
argued that aging might enhance the rate of evolution 
by decreasing generation overturn times |24) , reduce epi- 



demic outbreaks by diluting the population [35J or de- 
crease kin competition between parent and progeny |26j . 

While the mechanistic theories are firmly supported 
by experiment and continue to provide insight, the evo- 
lutionary origins of aging remain a mystery. The neutral 
(i.e. non-selective, non-directional) MA theory has two 
predictions: a monotonic increase in the mortality rate 
with age, and an increased variation (spread) in mortal- 
ity rate among different polymorphisms with age, both 
of which disagrees with observation [I7H19) . The non- 
neutral (i.e. selective, directional) AP theory predicts 
that every aging gene comes with an early-life enhance- 
ment of fecundity. While some such genes have been 
found, others contradict this prediction pZ0H23] . 

In contrast, the ideas of GS remain untested. Even if 
true, these theories do not explain why the mortality rate 
should increase in time, let alone why so many different 
organisms should exhibit a similar functional form for 
/i(i). 

Here we approach the question of aging from an evolu- 
tionary perspective, combining ideas from network con- 
nectivity [30H33] . engineering reliability analysis as ap- 
plied to aging [Ml [35], and the theory of constructive 
neutral evolution [27H29] . Much like the historical devel- 
opment of an engineered technological device, the com- 
plexity of life appears to have irreversibly increased as 
a large number of individual sub-units become linked 
through specialized interdependence. Our simulations 
and analysis allows us to demonstrate that component 
dependencies established over the course of neutral and 
non-neutral evolution eventually leads to species with ag- 
ing curves that are consistent with empirical data. 



II. ALGORITHM FOR NETWORK EVOLUTION 

To study the dynamics of aging quantitatively, we start 
with a simple view of an organism as a set of nodes 
with dependencies characterized by directed edges be- 
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FIG. 1. Fraction <j>(i) of functional components at 
time t as a function of dependency network size and 
topology, 100 Runs. Dependency structures grown via 
non-neutral (left column) and neutral (right column) evolu- 
tionary schemes yield scale free -P(fc) ~ 1/k 3 and random 
P(k) ~ e~ h / K degree distributions respectively. A represen- 
tative network for each scheme is shown on top of the respec- 
tive columns; preferential attachment produces high-degree 
hubs (large, black nodes) that are absent using random at- 
tachment. We plot 100 runs of <j>(t) for each network size and 
topology and show their respective lifetime r distributions 
/(t) in the inset. Increasing the network size from N = 2500 
(purple and pink) to N = 10 6 (blue and red) sharpens fir). 
The dashed black lines mark analytical predictions (cf. the- 
ory section) for initial slope po (which is 1.80 for SFN and 
1.75 for RN) and critical fraction <f>c (which is 0.6 for SFN 
and 0.5 for RN) which agree well with simulation results. For 
both plots {7o,7i,d} = {0.0025,0,0}. Note the remarkable 
similarity of different network topologies. 



tween them. Each node may be thought as genes in a 
regulatory network, or the differentiated cells or tissues 
in a multicellular organism with specific functions. A 
directed edge from node A to node B indicates that A 
provides something to B such as energy, crucial enzymes 
or mechanical support, so that the function of B relies on 
the function of A. In this scenario, the evolution of the 
network is represented by random addition of new nodes 
and edges, leading to a change in the dependency struc- 
ture of the network. This then changes the susceptibility 
of the network to further changes, including its longevity. 

For a given network, we assume that its dynamics are 
governed by four parameters, of which three are subject 
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FIG. 2. Empirical mortality data (black) fit 
to present theory (red) The mortality curves of 
(left to right) C. elegans, Drosophila, Medflies, Bee- 
tles, Mice, Himalayan Goats from [TH3]. Fit parame- 
ters are {JV.To.-yi.d} = ' t 700 . o.oi, 0.04, 0.025}, {800,0.0035,0.002,0.035}, 

{700, 0.01, 0.035, 0.015}, {500, 0.008, 0.065, 0.057}, {6000, 0.002, 0.014, 0.08}, 

{2100,0.038,0.29,0.096}. The horizontal axis (as well as the units 
of 1/70,1) is time in units of days for C. elegans, Drosophila, 
medflies and beetles; days/10 for mice and years for Tahr. To 
demonstrate that aging is a manifestation of component in- 
terdependence we also fit our model to two kinds of cars (cf. 
supplemental information and Fig.S.l) 



to direct experimental control: The failure rate 70 <C 1 
and repair rate 71 <C 1 of individual nodes, the initial 
fraction of damaged nodes d <C 1, and total number of 
nodes N 3> 1. 70 is controlled by biomolecular processes 
subject of mechanical theories of aging (e.g. oxidative 
stress, radiation damage); the damage due to pre or post 
natal stress is contained in d. 71 depends on the activity 
or inactivity of genes or their regulators that may be 
relevant for repair and replacement of cells. It seems 
difficult to determine or modulate N, though it should 
roughly correlate with the "complexity" of an organism. 

The structure of networks is expected to influence the 
dynamics of any processes on it, including aging. How- 
ever we surprisingly observe that different evolutionary 
processes for network growth lead to similar average lifes- 
pans and mortality curves. At either extreme of network 
structure, we use two strategies to grow them via what 
may be termed neutral and non-neutral evolution. In 
either case, a new node is introduced at every "evolu- 
tionary time step" such that it depends on one random 
existing node, and one random existing node depends on 
it. In the constructive neutral process, the new node is 
equally likely to depend on and provide to any of the ex- 
isting nodes. This yields a random dependency network 
(RN). In the constructive non-neutral process the proba- 
bility that a new component will provide or receive from 
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another component with degree k is taken proportional 
to k. This yields a scale-free dependency network (SFN) 
[55] , In the first scheme there is no difference in fitness 
based on whom a node relies on; whereas in the second 
scheme, relying on a previously unrelicd component has 
a larger fitness cost, making such connections unlikely. 
Each network topology is diagrammed schematically in 
Fig. [1]. 

These two schemes are not altogether arbitrary. They 
can be justified using biological arguments and are bol- 
stered by empirical evidence. SFN structures are known 
to accurately describe the organization of the metabolic 
networks of many different species [37] • A simple evo- 
lutionary mechanism [29] can give rise to a RN struc- 
ture and account for the complexity of various biological 
structures via an entropic argument: Since the number 
of configurations in which components randomly depend 
on each other is overwhelmingly more numerous than 
the configurations in which all elements arc fully inde- 
pendent, an initially- fully- independent collection of com- 
ponents will, statistically speaking, inevitably move to- 
wards random interdependence, provided that the fitness 
cost of introducing a dependency is negligible. 

We assume that the aging of complex dependency net- 
works are governed by the following three rules: (1) Ev- 
ery component in the organism must depend on at least 
one other node, and at least one other node must depend 
on it (i.e. all parts of the organism must be fully con- 
nected). (2) With certain fixed small probabilities the 
components can break (stop functioning) or be repaired 
(start functioning). (3) A node stops functioning if the 
majority of those on which it depends (providers) stop 
functioning, and cannot be repaired without a majority 
of its providers functioning. 

Our simulations are based on the three rules listed 
above implemented as follows. 

1. Create a network model of an organism 

(a) Begin with a single node, and i = 1. 

(b) Introduce a new [i + l) th node and make it 
depend on any one of the pre-existing nodes 
j < i with probability P(kj), where kj is the 
degree of node j. For the neutral scheme 
P(kj) is taken to be uniform and indepen- 
dent of kj , whereas for the non- neutral scheme 
P{kj) is taken proportional to kj. 

(c) Make any existing node j depend on the (i + 
l) th node with probability P(kj) 

(d) Increment i and repeat step b and c for N — 1 
steps. 

2. Age the resulting network model of an organism 

(a) Define the organism functional vector ip(t) = 
{xi(t), a:2(t), . . . , xjv(t)} where every compo- 
nent Xi can take either one of the values 1 
(functional) or (non functional). The vi- 
tality of the organism is defined as <j)(t) — 



YnXi(t)/N. Assign a value of to a fraction d 
of randomly selected nodes and 1 to the rest, 
corresponding to the initial damage of func- 
tionality in an organism. 

(b) For all i, update Xi = 1 to Xj = with proba- 
bility 70, flip Xi — to Xi — 1 with probability 
71, and do nothing with probability 1—70—71- 

(c) Break a node if the majority of nodes on which 
it depends are broken. Recursively repeat un- 
til no additional node breaks. 

(d) Set ip(t + 1) to the outcome of step (c). 

(e) Increment t and repeat (b-d) until all nodes 
are broken (i.e. Yli x i(t) = 0)- 

In order to study the network mortalities, we must 
define a time of "death" r, which could be chosen in 
multiple ways. We define a threshold 77 = 4>(t) = 1% 
below which an organism is defined dead [32] . The effect 
of the threshold value in the vitality <f>(t) that we use to 
define death does not affect our simulations for large N. 

To establish the statistical properties of mortality in 
our network, we generate an ensemble of networks (or- 
ganisms) and age them according to the above rules. In 
addition to tracking the vitality of the network character- 
ized by </>(i), we also determine the fraction of networks 
(not components) s(t) that remain alive at time t so that 
the time dependent mortality rate is 

fi(t) = -[s(t + l)-s(t)]/s(t) 

We also track the degree of interdependence within a 
given network (organism), characterized by the ratio 

A[0(*)]=log[0(*)]/]og[0o(i)], 

where 0o = exp{(— 70 + 71 )i} is the expectation value of 
the vitality of an identical size network with all depen- 
dency edges removed. Thus A quantifies how much more 
often a network dies in comparison with one that has no 
interdependent components. 

Finally, in order to analyze the magnitudes of function- 
ality loss we consider the probability distribution S[A(j)] 
of event sizes Acj> (i.e. changes in <fi). Each event repre- 
sents an individual disease or recovery, the final one of 
which is death. 



III. RESULTS 

In our simulations a typical organism starts its life by 
a slow decay of <f> at a rate of (0)70, where the dimen- 
sionless number (a) — 1.80 for scale free networks and 
(a) = 1.75 for random networks. As an increasing num- 
ber of nodes die, the system approaches a critical vitality 
4>(t) = <j) c when all live nodes suddenly collapse. Typi- 
cal trajectories for both network topologies as well as the 
values of (a) and <j) c are shown in Fig[l]. 
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FIG. 3. Mortality Rate p(t) = —dts/s as a function of time. In each panel we test the effect of one single parameter, keeping 
the others constant. While random networks (open markers) seem to age slightly faster compared to a scale free networks (filled 
markers), network topology does not seem to have a significant effect on the qualitative features of p, in line with experimentally 
observed universality of mortality curves of different species. We average over 100 networks with 1000 simulations each; thus 
the lowest probability event we can resolve is of the order ~ 10 -5 , and fluctuations on that order is likely noise, (a) A higher 
damage probability 70 shifts the lifespan distribution and mortality curve to the right (b) Repair rate changes the plateau 
value fio (c) Increasing N increases the slope of fi in the aging (Gompertz) regime. Only for large, complex networks do we 
find /i varying significantly with t; simple organisms do not age. (d) The initial damage causes a high infant mortality, but the 
damage is efficiently repaired soon after birth 



We observe that even if 71 is set equal to 70 the sys- 
tem decays steadily despite the seeming reversibility in 
dynamics. This is because while any live node can break, 
not all the dead nodes will have the sufficient number of 
live providers to sustain a repair. 

A number of qualitative features of the mortality 
curves predicted by our model are independent of the 
range of parameters we explore and the network struc- 
tures we chose. In particular, we see that the following 
features are generic (1) There is a slightly higher infant 
death rate proportional to initial damage d followed by 
a sudden drop in death rate in early childhood (2) There 
is an exponential increase in the mortality following this 
initial blip, consistent with the Gompertz Law and (3) 
There is a final plateau in late-life mortality. In our sim- 
ulations aging decelerates in late life and in any case does 
not increase beyond a critical po, the value of which is 
modulated by 71 (Fig[3b]). This is particularly clear in 
Fig [2], where we compare the mortality curves generated 
by our digital populations to that of a variety of organ- 
isms, C. elegans, Drosophila, Medflies, Beetles, Mice, Hi- 
malayan Goats (Tahr), using data compiled from [THS], 
and see reasonable agreement between simulation and 
data. 

Since complex interdependence in networks is not ex- 
clusive to living organisms, we also fit the empirical aging 



curves to our model, of two kind of automobiles, the 1980 
Toyota and 1980 Chevrolet (see Supplementary Informa- 
tion, FigS.l) as gathered from [3] . 

An empirical analysis of historical human mortality 
data has shown that plotting the mortality rate p(ti) 
at some age t\ against that at infancy p(to) gives a "uni- 
versal" curve that nearly overlaps for many different so- 
cieties and historical periods [38 . In other words, for an 
arbitrary collection of people (presumably with different 
damage /repair rates and initial conditions), the mortal- 
ity rate at any two ages are correlated, with correlation 
coefficient p ~ 1. In our simulations we vary 7o,7i,c2 
while keeping N constant and determine correlation co- 
efficients between mortality rates from ages tO to tO + 5 
and from t\ to t x + 5. For {to, t x } = {0, 20}, {10, 20} and 
{15,20} is correlated by p ~ 0.8,0.98 and 0.99. These 
correlation values are nearly identical for scale free and 
random networks. For details, see supplementary infor- 
mation and fig.S.6. 

The effects of the system parameters on the mortal- 
ity rate for both scale free and random networks can 
be summarized as follows: Increasing 70 shifts p(t) left 
(Fig[3a]); increasing 71 decreases the value of the value 
of po = — * 00 ) at the late life plateau (Fig[3b]); 
increasing N increases the slope of p in the aging (Gom- 
pertz) regime (Fig [3c]); in other words, larger systems 
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FIG. 4. Lifetime, interdependence and event distribu- 
tion, (a) Average lifetime (r) versus damage rate 70 (blue circles) and 
repair rate 71 (red squares) for scale-free networks (filled markers) and 
random networks (open markers) with TV — 2500. Dashed lines mark 
r — O.2O/70 for both random and scale free networks. 70 — 0.0065 
is kept constant for both red curves. There seems to be a critical 
repair rate 7* for which expected lifespan diverges. Note the remark- 
able independence of the curves with respect to network structure, (b) 
X(t) — log[0(£)]/7o* as a function of t (left) and cp (right) for 100 
trajectories on the scale free (blue) and random (red) networks with 
N = 10 6 pictured in Fig[l] (71 — 0). The interdependence parame- 
ter A varies strongly with both t and and roughly doubles as more 
damage is accumulated, until the sudden collapse depicted in Fig[l] 
leads to a diverging interdependence. Scale free networks show a rel- 
atively large variation in A for short times, but rapidly converge on a 
monotonic increase as the network accumulates damage. Interestingly, 
scale free networks begin with a larger value of A, but random networks 
become more interdependent rapidly. The interdependence of a set of 
disconnected nodes (completely independent) have Ao — 1 shown in the 
dashed line, (c) Probability S that drops by At/> before the largest 
drop happens (blue) of the largest drop itself (purple) and after the 
largest drop happens (yellow). Scale free (solid) and random (dashed) 
networks show a remarkable similarity. The largest drop distribution 
of both random and scale free networks obey a power law with ex- 
ponent — 2.7, (black line marks slope). Note that "disease" (blue) is 
qualitatively different from "death" (purple). Simulation parameters 
arc {N, 7o,7l, d} = {2500,0.0025,0,0}. 



age rapidly and suddenly, while small systems with few 
components are virtually non-aging. Increasing the ini- 
tial damage d simply elevates the initial (infant) mortal- 
ity rate (Fig[3d]). Our simulations yield a negative cor- 
relation between damage and increase in mortality rate; 
the high initial damage populations age slower than the 
low initial damage populations to eventually converge to 
the same /j,q consistent with Strehler-Mildvan correlation 
law @T]. 

The qualitative dependence of average lifetime on dam- 
age and repair rate is as intuitively expected (Fig [4a]). 
Quantitatively, when 71 = the average lifespan per- 
fectly fits the curve (r) = /3/7o, with the same value 
of j3 = 0.2 fitting both scale free and random networks. 
Curiously, lifespan is very weakly dependent on the re- 
pair rate 71 for smaller values until it rapidly diverges 
at a critical value of 7*. This makes one wonder why 
most species are not immortal, but Fig[4a] suggests why: 
For small 71, r(7i) is nearly constant, (\/T)dr/d~f\ <C 1. 
However, if the cost of repairing nodes increases with the 
probability of repair 71 , increasing 71 becomes evolution- 
arily nonviable due to the weak dependence of t over a 
wide range of 71. That being said, apart from this high 
evolutionary cost (and the slim chance of a massive sta- 
tistical fluctuation of the order 7^ which obviously van- 
ishes for N — > 00) there is no mechanism in our theory 
that prevents a high-repair species from being immortal 
in the thermodynamic limit. 

The dependency coefficient A[0(i)] gradually increases 
as our simulated organisms grow older, and diverges just 
before death (Fig[4b]). The functional form of A[t/>(i)] is 
qualitatively similar for SFN and RN. 

We intuitively expect the onset of "death" to differ 
drastically from the early aging process (referred to as 
"disease" ) . This difference will be reflected in the distri- 
bution of the number of living nodes that die in a par- 
ticular time step (Ac/) = cf>(t) — <j)(t — 1)). In order to tell 
whether death is just "the last disease" or a qualitatively 
different phenomenon, we analyze the distribution S(A<f>) 
of event sizes before, after and of the largest drop and 
notice that the latter is qualitatively as well as quantita- 
tively very different from the former two. Death and dis- 
ease occupy an entirely different region of the event spec- 
trum (Fig[4c]). What is even more remarkable is that the 
disease distribution for both evolutionary schemes RN 
and SFN are quantitatively similar over a wide range of 
A<p, and obey a power law S(A(f>) ~ l/A(f> 2 - 7 . We do not 
have an explanation for this striking similarity, nor the 
value of the critical exponent. 

We cautiously note that while a set of four parameters 
uniquely determines /i(i), the converse is not true. In 
particular the characteristic features of /1 (as defined by 
the initial slope, the plateau value /in,, average lifespan r 
and crossover time from an aging to non-aging regime) 
can be kept "similar" by holding 70 constant while in- 
creasing N, 71 and d simultaneously (see Supplementary 
Information and Fig.S.4). Of course, there is no a pri- 
ori reason why two species with different attributes must 
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necessarily have different aging curves. 



IV. THEORY 

We now aim to obtain the values of initial decay rates 
(a)70i the critical vitality <p c and understand why depen- 
dency networks collapse suddenly. On the way, we also 
hope to understand why these quantities are so similar 
for both scale free and random networks, and determine 
the origin of the Gompertz-like law. 

When the system is far from collapse, the probability 
that two providers of a single node dying at once 0[^q] 
is negligible compared to that of a single provider dying 
O[7o]. Then the total probability po that a node dies 
is 70 plus the probability that the last vital provider of 
a node dies. If m(<p) is the probability that a node is 
left with one last vital provider, we can self-consistently 
evaluate po 



Pa = 7o + ™(</>)Po(l - 7o)- 



(1) 



In a single step associated with the aging of the network, 
the probability that a node is repaired is p\ — h((j))^\, 
where h(<fi) is the probability that a node has at least 
the minimum number of providers required to function. 
Then, the change in the fraction of nodes that are alive, 
is given by 



A(j> = pot 



700 

l-m(0)(l- 7o ) 



+ 7iM0)(l-0) (2) 



where we have used the expression for p as obtained 
from (JlJ. From ([2]), we see the origin of the catastrophic 
(and universal) nature of death. For any arbitrary "fully 
connected" network and monotonically decreasing 4>{t), 
the vital fraction m(<p) must always start from a finite 
value in the domain [0, 1] and increases towards unity as 
<t> decreases, inevitably to cause the first term to dominate 
the second (70 -C 1), and thus leading to a sudden drop 
in the expected vitality. This is true in general, although 
the detailed form of the evolution of <fi depends on the 
fraction m(4>) of vital providers and repairable fraction 
h((j)) that will vary for different network structures. 

Equation ([2| also indicates an asymptote in longevity 
for large repair rates, as seen in Fig [4a]. If we set A<fi = 0, 
we find that the system lives indefinitely when the repair 
rate is set to 71 = 7* where 



7i 



7o</>* 



ft(0*)(l-0*)[l-m(0*)(l-7b)] 



(3) 



for any given </>* £ [<p c , 1], i.e. for a vitality larger than 
the critical vitality <f) c . In this case the system damage 
increases while the vitality decreases until it reaches cf> = 
cf>* , but maintains that damage forever. Of course, ^ 
governs the expectation value of <f>(t) which is the actual 
value only in the thermodynamic limit N — > 00. Thus, 



only a system of infinite size that satisfies ^ can live 
indefinitely, since a finite system will die (at least) with 
probability -y^^'*J due to statistical fluctuations. 

In general, it is a non-trivial task to obtain the exact 
forms of m((j)) and h((f>) and thence the average lifetime, 
the critical damage fraction etc. However we can obtain 
the initial slope (0)70 = pa\t=o with which the vitality 
decreases, and the critical vitality C at which the whole 
system collapses (see dashed lines in Fig[l]) for the case 
71 = 0: Let the probability that a node with k providers 
die be a^ k \ Then we can recursively obtain cr*- 1 -* in terms 
of the others HO] , 



a« = 7 o + P(l,l) ( T (1) 



P(l,2)tr< 2 > +P(l,3)cr 



(3) 



(4) 



where P(l,«) is the probability that the provider of a 
degree- 1 node has i providers. The first term corresponds 
to the probability that a node dies independent of its 
connectivity. 

Since we neglect probabilities of order 0[7q], initially 
only degree-1 nodes can be killed by the death of their 
providers. Thus substituting = 70 for all k apart 
from k = 1, and using $~1-P(l,i) = P(l) we can obtain 
from Q the initial probability that a degree-1 node dies, 



r(i) - 



(2-P(l,l)) 70 
l-P(l.l) 



(5) 



to find the expectation value of the initial slope we must 
average over the damage rate of all degrees, including 
o-W = 70 for k > 1 



{a) la 



(it) 



7o 



t=a 



P(l) 



1-P(1,1) 



(6) 



Upon substituting the numerical values of P(l) and 
P(l, 1) for the networks we evolved, we obtain (a) = 1.75 
for the neutral scheme (RN) and (a) = 1.80 for non- 
neutral (SFN) scheme, i.e. only a ~ 2.8% difference be- 
tween two topologies. These initial slopes are consistent 
with our aging simulations (Fig[l]). 

To estimate the critical vitality <f> c , we start by asking 
what the smallest value of 70 must be in order to kill any 
network in just one step. If there exist a critical vitality 
4> Cl then a death rate of 70 = 1— <f> c would kill the network 
in one step. Making this substitution and letting o"W — > 1 
for all i in Q yields a simple but interesting result, 



0c = P(l). 



(7) 



For the networks we evolved, P(l) is 0.5 for the neutral 
scheme and 0.6 for the non-neutral scheme. These values 
are consistent with the critical vitalities we observed in 
our aging simulations (Fig[l]). 

Having estimated the average damage rate and the 
critical vitality of the network, we now consider the na- 
ture of lifespan distributions. The probability of network 
survival s(t) is equal to the probability of 8<f> being not 
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greater than (j> — <fi c . That is, s(t) = 1 — Prob(A</> > 
4> — <fi c ). Since each node dies with probability po-, 

n\ 1 ( N< t ) \ N(4>-4> c )+k n \N<f> c — 

i(t)=1_ SU-w+*J ft ( po) 

(8) 

Using the aging rate of the network established earlier, 
i.e. 1.757o < Po < 1 regardless of the network topology 
or size, i.e. po is independent of both these parameters. 
Thus, in the limit N — > 00 the probability of death 1— s(t) 
simply becomes a unit step function. On the other hand 
for finite N, the step function softens, and we empirically 
interpret the rapid transition from s = 1 to s = as 
aging. We see that by passing from finite size to infinite 
size, we also pass from the stochastic to deterministic 
case, and from gradual aging to "instant aging". 

Finally, we consider the sharpness of the transition 
from s = 1 to s — to see if there is any relation between 
our results and the classical Gompertz law for mortality. 
Evaluating ^ exactly is nontrivial, since, strictly speak- 
ing, all quantities in ([sj) except N are random variables. 
However we can substitute the average (early-time) value 
of <p(t) = e~ Pot in this is done at the cost of narrowing 
the longevity distributions but hopefully not their func- 
tional forms. The mortality rate /j,(t) — —dt-s/s obtained 
this way is plotted in supplemental Fig.S.2 and agrees 
well with the empirical Gompertz law, which states that 
log/x(t) increases linearly in time. 

V. DISCUSSION 

We have built on a rough similarity between large net- 
works and complex organisms (and machines) to create 
a minimal model of aging. It is thus important to be 
self-critical here by comparing our study both with real- 
ity and with previous attempts. Any theory that aims to 
account for a phenomenon as universal as aging spanning 
both animate and inanimate objects needs to be robust, 
i.e. it should not have such strong assumptions and re- 
sults that strongly depend on system details. The weak 
sensitivity of our outcomes to the details of biologically 
justifiable dependency network structures (Fig[l,2,3,4]) 
seems to satisfy this requirement. However, as already 
pointed out, we have four parameters in our model - 
network size, rate of damage and repair, and the ini- 
tial damage, in addition to the fraction of dependency 
that determines if a node dies or not (1/2 in our case). 
It seems unlikely that any model could be simpler, and 
while this is sufficient to fit a variety of mortality curves 
in organisms and machines, the set of these parameters is 
not unique for an organism. Clearly, therefore, we need 
further constraints to make our model even more robust, 
and some of these may come from mechanistic limits on 
the parameter values. 

Comparing our analysis with earlier models, we note 
that a classical benchmark is provided by the application 



of reliability analysis to biological aging [Ml |3"5] , which 
investigates the robustness of systems in which cells are 
connected in a parallel-series circuit (i.e. an organ lives 
as long as one of its cells remains, and the organism dies 
as soon as one of its organs fails) . While this picture can 
successfully explain some qualitative features of mortal- 
ity curves, a few puzzles remain. First, it cannot account 
for the higher infant mortality followed by a drop (in con- 
trast see Fig[3d] and second row of Fig. S. 6) despite intro- 
ducing an initial damage. Secondly, the plateau value of 
/io = hnit^co jit(t) implied by the model of [M] regardless 
of the number of cells per organ, is equal to 1 — (1 — 70 ) m , 
where m is the total number of organs. This implies that 
any link to actual experimental values of fj,o for an organ- 
ism such as C. elegans, with say, m ~ 10 organs, requires 
about ten percent of all cells die per day, an unreason- 
ably large fraction that contradicts experiment. The final 
puzzle in [34, 35] is, the predicted functional form of /i(t) 
depends very sensitively on a specific assumption of how 
the reliability circuit is built. Specifically, the functional 
form of the probability distribution describing the num- 
ber of "redundant cells per organ" which by itself seems 
to have no a priori or experimental justification. In con- 
trast, our analysis of random and scale free networks is 
able to explain the initial increase in mortality, is consis- 
tent with "low" rates of cell deaths seen in experiments, 
and is surprisingly independent of connectivity distribu- 
tions or individual realizations of networks. 

Our analysis also differs significantly from earlier the- 
oretical investigations of network failure 30-33 in which 
nodes are simply removed one by one (systematically or 
randomly) until networks get fragmented. In these ap- 
proaches the nodes do not influence the performance of 
one other, and are not allowed to be repaired. Curiously, 
the lack of interactions in these models lead to fundamen- 
tally different fragmentation dynamics in scale free and 
random networks. In contrast the survival curves we ob- 
serve are remarkably independent of the network topol- 
ogy (Fig[3,4]). The strong interactions between compo- 
nents may be the reason behind the strong similarity of 
mortality curves among so many organisms; a conclusion 
that bears similarity to that of another strong-interaction 
model that has been proposed to explain electrical gird 
failures [33] . Much still remains to be done in under- 
standing how the form of these interactions leads to dif- 
ferences or similarities in the dynamics. 

Our study has focused on the dynamics of networks 
that age as a consequence of interdependency, and thus 
leads naturally to the question of how this might be con- 
trolled. Since the repair rate, and perhaps the damage 
rate (to a lesser extent) are experimentally controllable, 
one might ask if it is possible to vary their temporal char- 
acter while keeping their average constant. Arc there op- 
timal strategies for repair - either in the time domain or 
in space (i.e. looking at nodes with varying connectiv- 
ity)? For example, Fig [4a] shows that if the system is 
repaired uniformly, the degree of repair does not make a 
significant difference for 71 < 7* . It would be very useful 



to know if and how a (temporal or spatial) non-uniform 
repair strategy improves lifespan. In networks that are 
dynamically heterogeneous, we may ask what would be 
the consequences of differential damage and repair in a 
network with highly variable turnover - e.g. a network 
with tissues like the skin or gut that have high damage 
and repair rates, and the brain which has a low damage 
and repair rate? At the level of ecology and colonies, we 
might ask how does the aging dynamics of a dependency 
network change when a system consists of parts with ag- 
ing rates comparable to that of the whole system? We 
hope that our minimal model may be used to study some 
of these questions. 
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the random variables <p(t) and po by their average value 
we sharpen the lifetime distributions and hence steepen 
the mortality curve; however Gompertz's Law is still 
recovered for short times (Fig.S.2). 

3. We have defined death as the time r at which <f> 
reaches a threshold value 1%, which may seem arbitrary. 
In order to determine how sensitive our results are to 
the choice of threshold rj, we have analyzed the lifetimes 
of both network topologies as a function r\. For a 
(small) network of /V = 2500 Fig.S.3 shows that the 
value of r\ changes the lifetime less than 1% for a SFN 
and less than 12 for a RN. We find that r\ dependence 
rapidly vanishes for N > O[10 3 ] for both network types 
and conclude that the precise value of r\ is not important. 

4. To demonstrate that our model parameters have 
individual predictive power we fit the empirical data of 
a mutant and wildtype of a fixed organism C.Elegans by 
fixing iV,7o and d constant, but only varying the repair 
rate 71 (Fig S.4). The difference in fit parameters as 
noted in the legend of Fig[2] and Fig. S.4 stem from a 
discrepancy in the data presented in [3] and p] for the 
same organism. 



Appendix: Supplemental Information 

Here we discuss in more detail various technical points 
surrounding theory, simulations, and fits. 

1. Although we have discussed aging in an evolution- 
ary setting, specialized interdependence of components 
is not exclusive to living organisms. To demonstrate this 
point we fit our model to the empiric mortality curves 
of 1980 Toyota and 1980 Chevrolet obtained from [3]. 
We were able to obtain a reasonable agreement using 
identical values of N and d but slightly varying damage 
and repair rates for the two cars (Fig S.l). 




t t 



FIG. S.l. Empiric mortality rates of cars (black) fit to 
present theory (red). The data for 1980 Toyota (left) and 
1980 Chevrolet (right) from [3] is fitted with {N, 70, 71, d}= 
{200, 0.023, 0.023, 0} and {200, 0.02, 0.02, 0} respectively. The 
horizontal axis denotes years. 

2. In order to recover the Gompertz Law analytically, 
we substitute 4>(t) ss e ~ a7ot for t <C t in (8) and plot 
/i(<) = —(ds(t)/dt)/s(t) in Fig.S.2. By approximating 
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FIG. S.2. Mortality vs time obtained from (8) for cf> <C 4> c , 
70 = 0.0025 and N = 50 (black) 100 (red) and 200 (blue) is 
in good qualitative agreement with our simulations and the 
empiric "Gompertz law" which states that log fi(t) increases 
linearly with t early in life. 



5. To determine whether model parame- 

ters {N, 70, 71, d} can be determined uniquely 
given experimental /i(i) data we performed sim- 
ulations sweeping the parameter space N € 
{50, 100, 250, 500, 700, 800, 900, 1000, 1500, 2000}, 
7o G {2,4,6,8,10,15,20,25,30,35} x 10" 3 ,7i e 
{1,2,..., 10} x 7o,d € {0,0.5,...., 10} (percent) and 
checked if non-neighboring parameters give more similar 
n(t) curves than neighboring ones (cf. below for details). 
For each set of parameters we generated 12 networks 
upon which 3000 simulations were performed, provid- 
ing a reasonable level of confidence in the statistical 
accuracy of the simulations. To quantitatively compare 
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FIG. S.3. Sensitivity of outcomes to the definition of 
death, (a) The percent lifetime difference between choices 
r\ = 1% and 50% for scale free (filled circles) and random 
(empty circles) dependence networks, with varying 70 and 
71 = 0. There is a fairly constant 5% difference independent 
of 70 for random networks, and below 1% for scale free net- 
works, (b) The percent lifetime difference between choices 
■q = 1% and 50% for varying 71, with 70 = 0.00625. Scale free 
networks continue to have below 1% difference between the 
two threshold choices, while the variation is more significant 
for random networks with large repair rate 71. The relative 
difference decreases if the two thresholds r\ are closer to one 
another. For both graphs N — 2500, d = 0. 



the simulation results, (i(t) is broken into four averaged 
characteristics: The initial slope, the saturation point, 
the crossover time between the initial growth and 
saturation, and the observed lifetime (see Fig. S. 5a). The 
threshold for similarity of the curve characteristics is 
determined by averaging over the differences in nearest 
neighbors in the (TV, 70,71, d) parameter space (i.e. the 
5Zfc=i \ T ref ~ T s k |, with the parameters in the simulation 
Sk being different from the reference simulation in only 
one position, and a nearest neighbor). Thus, simulation 
outcomes are considered "similar" to a reference if they 
are not nearest neighbors (in parameter space) with the 
reference, and if the differences in all four characteristics 
simultaneously fall within the threshold variation. 
Fig. S. 5b shows one such overlapping curve with d = 4.5, 
N = 700 and 71 = 570 (compared to d = 0%, TV = 250 
and 71 = for the reference simulation in black). The 
inset of Fig. S. 5b shows all 214 simulation parameters 
that yields mortality curves considered "similar" to 
the reference curve (about 0.9% of all of the simulated 
parameters), and shows that a rather wide range of 
parameters may give qualitatively similar behavior in 
fi and t (with the latter not shown). It is interesting 
to note that 70 is the same for the reference curve 
and all overlapping curves, indicating that an empiri- 
cally observed death rate 70 may be uniquely determined. 

6. It has been empirically observed that the value of 
n{t{j) at age to correlates very strongly with /it(ti) at age 
ti for a wide variety of societies and historical periods 
[55] , To test whether our model yields this empiric 
"universality", we plot fi(to) versus n{t\) for a varying 
range of 70,71, d and network types. While the trend 
is less impressive for very widely separated t a ,tl, our 
model does yield a correlation between rait pairs for a 
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FIG. S.4. Mortality Rate of Mutant and Wild- Type 
Nematodes. We fit the data from Fig. 3E of [3] using our 
model. The data includes a wild type C. Elegans, which is 
well fit by the parameters TV = 500, 70 = 0.018, 71 = 270, 
and d = 0. The mortality curve for the Age-1 mutant of the 
worm is well fit using the same values of TV, 70, and d, but 
with 71 = IO70 significantly increased. 
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FIG. S.5. Uniqueness of fj,(t) as determined by initial 
slope, saturation value crossover time and average 
lifetime (a) A simulation is quantified in terms of four pa- 
rameters: The initial slope, the final saturation value, the 
crossover time, and the lifetime. Shown is a reference simula- 
tion with TV = 250, 70 = 0.002, 71 = 0, and d = 0% (b) The 
non-uniqueness of fi(t) as the parameters are varied. In the 
main panel and the inset, the black points correspond to the 
reference simulation in (a). The three red line in the main 
panel has TV = 700, 70 = 0.002, 71 = 570, and d = 4.5%. The 
average lifetime for the red curve is r = 133, within 4.3% of 
the lifetime of the reference curve. There is moderate vari- 
ation between the curves for small and large t, but it would 
be difficult to unambiguously differentiate between the two 
sets of parameters when fitting experimental data. The inset 
shows the same reference simulation (black points), along with 
all 214 sets of simulated parameters that satisfy the threshold 
criterion. Each blue line has 250 < TV < 2000, < 71/70 < 9, 
and < d < 8.5, all with 70 = 0.002. 
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FIG. S.6. Testing for universality in mortality patterns. Dependence of old age mortality fj,(5) = Prob(20 < r < 25) on 
younger age mortality fi(i) — Prob(5i — 5 < r < 5i) for i — 1 (top left), 3 (top center), 4 (top right) for a range of network 
types (RN blue, SFN red) and system parameters, 70 = {0.0026,0.0027,0.0028,0.0029}, 71 = {0.0026,0.0027,0.0028,0.0029}, 
d — {12.1, 12.2, 12.3, 12.4} (r is the time of death). The mortality curves generated by this range of parameters is displayed in 
the bottom row, for RN (left) and SFN (right). While the universal (yet species-specific) trend observed in [38] between fi(i) 
and fj,(i + j) (for fixed j) is qualitatively present in our model, the trend vanishes for large enough j. 



range of network parameters and network types (fig.S.6). 
The correlation coefficients between mortality rates from 
ages tO to tO + 5 and £i to h + 5 is p ~ 0.8, 0.98 and 0.99 
for {tQ,h} = {0,20}, {10, 20} and {15,20} respectively. 



When determining these values the bin size was chosen 
as 5 instead of 1 in order reduce finite size effects (e.g. 
not many networks die exactly on step-1). 
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